/*
 * Copyright (c) 2019, Oracle and/or its affiliates. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *  http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

package org.eclipse.imagen.media.opimage;

import java.awt.Rectangle;
import java.awt.geom.Point2D;
import java.awt.image.DataBuffer;
import java.awt.image.Raster;
import java.awt.image.RenderedImage;
import java.awt.image.WritableRaster;
import java.util.Map;
import org.eclipse.imagen.BorderExtender;
import org.eclipse.imagen.GeometricOpImage;
import org.eclipse.imagen.ImageLayout;
import org.eclipse.imagen.Interpolation;
import org.eclipse.imagen.InterpolationBicubic;
import org.eclipse.imagen.InterpolationBicubic2;
import org.eclipse.imagen.InterpolationBilinear;
import org.eclipse.imagen.InterpolationNearest;
import org.eclipse.imagen.RasterAccessor;
import org.eclipse.imagen.RasterFormatTag;
import org.eclipse.imagen.media.util.ImageUtil;

/**
 * A class extending <code>GeometricOpImage</code> to subsample and antialias filter images. Image scaling operations
 * require rectilinear backwards mapping and padding by the resampling and filter dimensions.
 *
 * <p>When applying scale factors of scaleX, scaleY to a source image with width of src_width and height of src_height,
 * the resulting image is defined to have the following bounds: <code></pre>
 *       dst minX  = round(src minX  / scaleX) <br>
 *       dst minY  = round(src minY  / scaleY) <br>
 *       dst width  =  round(src width  / scaleX) <br>
 *       dst height =  round(src height / scaleY) <br>
 * </pre></code>
 *
 * <p>The applied filter is quadrant symmetric (typically antialias + resample). The filter is product-separable,
 * quadrant symmetric, and is defined by half of its span. For example, if the input filter, qsFilter, was of size 3, it
 * would have width and height 5 and have the symmetric form: qs[2] qs[1] qs[0] qs[1] qs[2] Because we have chosen to
 * keep the filters in compact form we need to keep track of parity.
 *
 * <p>A fully expanded 5 by 5 kernel has format (25 entries defined by only 3 entries): <code>
 *   <p align=center> qs[2]*qs[2]  qs[2]*qs[1]  qs[2]*qs[0]  qs[2]*qs[1]  qs[2]*qs[2] <br>
 *
 *                    qs[1]*qs[2]  qs[1]*qs[1]  qs[1]*qs[0]  qs[1]*qs[1]  qs[1]*qs[2] <br>
 *
 *                    qs[0]*qs[2]  qs[0]*qs[1]  qs[0]*qs[0]  qs[0]*qs[1]  qs[0]*qs[2] <br>
 *
 *                    qs[1]*qs[2]  qs[1]*qs[1]  qs[1]*qs[0]  qs[1]*qs[1]  qs[1]*qs[2] <br>
 *
 *                    qs[2]*qs[2]  qs[2]*qs[1]  qs[2]*qs[0]  qs[2]*qs[1]  qs[2]*qs[2]
 *   </p> </code>
 *
 * <p>Horizontal and vertical kernels representing convolved resample and qsFilter kernels are computed from the input
 * filter, the resample type, and because the downsample factors affect resample weights, the downsample scale factors.
 * If the scale factors are odd, then the resample kernel is unity. Parity is used to signify whether the symmetric
 * kernel has a double center (even parity) or a single center value (odd parity).
 *
 * <p>This operator is similar to the image scale operator. Important differences are described here. The coordinate
 * transformation differences between the FilteredDownsampleOpImage and the ScaleOpImage operators can be understood by
 * comparing their mapping equations directly.
 *
 * <p>For the scale operator, the destination (D) to source (S) mapping equations are given by <code>
 *   <p> xS = (xD - xTrans)/xScale <br>
 *       yS = (yD - yTrans)/yScale
 * </code>
 *
 * <p>The scale and translation terms are floating point values in D-frame pixel units. For scale this means that one S
 * pixel maps to xScale by yScale D-frame pixels. The translation vector, (xTrans, yTrans), is in D-frame pixel units.
 *
 * <p>The filtered downsample operator mapping equations are given by <code>
 *   <p> xS = xD*scaleX + (int)(hKernel.length/2) <br>
 *   yS = yD*scaleY + (int)(vKernel.length/2)
 * </code>
 *
 * <p>The mapping equations have the intended property that the convolution kernel overlays the upper left source pixels
 * for the upper left destination pixel.
 *
 * <p>The downsample terms are restricted to positive integral values. Geometrically, one D-frame pixel maps to scaleX
 * by scaleY S-frame pixels. The combination of downsampling and filtering has performance benefits over sequential
 * operator usage in part due to the symmetry constraints imposed by only allowing integer parameters for scaling and
 * only allowing separable symmetric filters. With odd scale factors, D-frame pixels map directly onto S-frame pixel
 * centers. With even scale factors, D-frame pixels map squarely between S-frame pixel centers. Below are examples of
 * even, odd, and combination cases.
 *
 * <p>s = S-frame pixel centers <br>
 * d = D-frame pixel centers mapped to S-frame <kbd>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>   d       d       d  </pre>
 *
 * <pre> s   s   s   s   s   s           s   d   s   s   d   s  </pre>
 *
 * <pre>  </pre>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>   d       d       d  </pre>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>  </pre>
 *
 * <pre> s   s   s   s   s   s           s   d   s   s   d   s  </pre>
 *
 * <pre>   d       d       d  </pre>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>  </pre>
 *
 * <pre> Even scaleX/Y factors            Odd scaleX/Y factors  </pre>
 *
 * <pre>   </pre>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>     d           d    </pre>
 *
 * <pre> s   s   s   s   s   s           s d s   s d s   s d s  </pre>
 *
 * <pre>   </pre>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>     d           d    </pre>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>   </pre>
 *
 * <pre> s   s   s   s   s   s           s d s   s d s   s d s  </pre>
 *
 * <pre>     d           d    </pre>
 *
 * <pre> s   s   s   s   s   s           s   s   s   s   s   s  </pre>
 *
 * <pre>   </pre>
 *
 * <pre>  Odd/even scaleX/Y factors      Even/odd scaleX/Y factors  </pre>
 *
 * <br>
 * </kbd>
 *
 * <p>The convolution kernel is restricted to have quadrant symmetry (qs). This type of symmetry is also product
 * separable. The qsFilter is specified by a floating array. If qsFilter[0], qsFilter[1], ... , qsFilter[qsFilter.len()
 * - 1] is the filter input, then the entire separable kernel is given by <br>
 * qsFilter[qsFilter.len() - 1], ... , qsFilter[0], ... , qsFilter[qsFilter.len() - 1] <br>
 *
 * <p>The restriction of integer parameter constraints allows full product separablity and symmetry when applying the
 * combined resample and filter convolution operations.
 *
 * <p>If Bilinear or Bicubic interpolation is specified, the source needs to be extended such that it has the extra
 * pixels needed to compute all the destination pixels. This extension is performed via the <code>BorderExtender</code>
 * class. The type of border extension can be specified as a <code>RenderingHint</code> to the <code>JAI.create</code>
 * method.
 *
 * <p>If no <code>BorderExtender</code> is specified, the source will not be extended. The scaled image size is still
 * calculated according to the formula specified above. However since there is not enough source to compute all the
 * destination pixels, only that subset of the destination image's pixels which can be computed, will be written in the
 * destination. The rest of the destination will be set to zeros.
 *
 * <p>The current implementation of this operator does not support <code>MultiPixelPackedSampleModel</code> source data.
 * The Rendered Image Factory for this operator, <code>FilteredSubsampleRIF</code>, will throw an <code>
 * IllegalArgumentException</code> for this type of input.
 *
 * @see GeometricOpImage
 */
public class FilteredSubsampleOpImage extends GeometricOpImage {

    /** The horizontal downsample factor. */
    protected int scaleX;

    /** The vertical downsample factor. */
    protected int scaleY;

    /** Horizontal filter parity. Rules: 0 => even, 1 => odd. See hKernel */
    protected int hParity;

    /** Vertical filter parity. Rules: 0 => even, 1 => odd. See vKernel */
    protected int vParity;

    /**
     * Compact form of combined resample and antialias filters used by computeRect method. hKernel is the horizontal
     * kernel.
     *
     * <p>The symmetric filter is applied even or odd depending on filter parity.
     *
     * <p>Expanded even kernel example (<code>hParity = 0</code>): <br>
     * <code>
     *      hKernel[2] hKernel[1] hKernel[0] hKernel[0] hKernel[1] hKernel[2]
     *  </code>
     *
     * <p>Expanded odd kernel example (<code>hParity = 1</code>): <br>
     * <code>
     *      hKernel[2] hKernel[1] hKernel[0] hKernel[1] hKernel[2]
     *  </code>
     */
    protected float[] hKernel;

    /**
     * Compact form of combined resample and antialias filters used by computeRect method. vKernel is the vertical
     * kernel.
     *
     * <p>The symmetric filter is applied even or odd depending on filter parity.
     *
     * <p>Expanded even kernel example (<code>vParity = 0</code>): <br>
     * <code>
     *      vKernel[2] vKernel[1] vKernel[0] vKernel[0] vKernel[1] vKernel[2]
     *  </code>
     *
     * <p>Expanded odd kernel example (<code>vParity = 1</code>): <br>
     * <code>
     *      vKernel[2] vKernel[1] vKernel[0] vKernel[1] vKernel[2]
     *  </code>
     */
    protected float[] vKernel;

    /**
     * <code>convolveFullKernels</code> -- convolve two kernels and return the result in a floating array.
     *
     * @param a floating kernel array.
     * @param b floating kernel array.
     * @return floating kernel array representing a*b (full convolution)
     */
    private static float[] convolveFullKernels(float[] a, float[] b) {
        int lenA = a.length;
        int lenB = b.length;
        float[] c = new float[lenA + lenB - 1];

        for (int k = 0; k < c.length; k++)
            for (int j = Math.max(0, k - lenB + 1); j <= Math.min(k, lenA - 1); j++) c[k] += a[j] * b[k - j];

        return c;
    } // convolveFullKernels

    /**
     * <code>convolveSymmetricKernels</code> uses a symmetric representation (partial kernels) of input and output
     * kernels. For example, with aParity 1 (odd) and bParity 0 (even) the passed kernels a and b would have form:
     * <code>
     * a:  a[lenA-1] ... a[1] a[0] a[1] ... a[lenA-1]
     * b:  b[lenB-1] ... b[1] b[0] b[0] b[1] ... b[lenB-1]
     * </code> (i.e., don't send symmetric parts but assumes parity controls filter lengths).
     *
     * <p>It is possible to do this convolution without resorting to full kernels but this is messy. @see
     * convolveFullKernels for details.
     *
     * <p>Further notes: 1. The return kernel, <code>c</code>, has parity <code>1 + aParity + bParity mod 2</code> 2.
     * The reason for setting up the kernels this way is to enforce symmetry constraints. (Design choice.)
     *
     * @param aParity int that is 0 or 1.
     * @param bParity int that is 0 or 1.
     * @param a floating partial kernel array.
     * @param b floating partial kernel array.
     * @return symmetric portion of floating array representing a*b (convolution).
     */
    private static float[] convolveSymmetricKernels(int aParity, int bParity, float[] a, float[] b) {
        int lenA = a.length;
        int lenB = b.length;
        int lenTmpA = 2 * lenA - aParity;
        int lenTmpB = 2 * lenB - bParity;
        int lenTmpC = lenTmpA + lenTmpB - 1;
        float[] tmpA = new float[lenTmpA];
        float[] tmpB = new float[lenTmpB];
        float[] tmpC;
        float[] c = new float[(lenTmpC + 1) / 2];

        // Construct "full" a
        for (int k = 0; k < lenTmpA; k++) tmpA[k] = a[Math.abs(k - lenA + (aParity - 1) * (k / lenA) + 1)];

        // Construct "full" b
        for (int k = 0; k < lenTmpB; k++) tmpB[k] = b[Math.abs(k - lenB + (bParity - 1) * (k / lenB) + 1)];

        // Convolve "full" a with "full" b to get a "full" tempC
        tmpC = convolveFullKernels(tmpA, tmpB);

        // Carve out and return the portion of c that holds
        int cParity = tmpC.length % 2;
        for (int k = 0; k < c.length; k++) c[k] = tmpC[lenTmpC - c.length - k - 1 + cParity];

        return c;
    } // convolveSymmetricKernels

    /**
     * <code>combineFilters</code> based on <code>resampleType</code> and input partial <code>qsFilter</code> (see above
     * for details on qsFilter format). Input <code>qsFilter</code> is restricted to have odd parity (<code>qsFilter[0]
     * </code> is at the center of the kernel).
     *
     * @param scaleFactor positive int representing the downsample factor.
     * @param resampleType int representing the interpolation type.
     * @param qsFilter floating partial kernel array (antialias filter).
     * @return floating partial kernel representing combined resample and qsFilter.
     */
    private static float[] combineFilters(int scaleFactor, int resampleType, float[] qsFilter) {

        // Odd scale factors imply no resample filter is required
        // return pointer to the qsFilter
        if ((scaleFactor % 2) == 1) return (float[]) qsFilter.clone();

        int qsParity = 1;
        int resampParity = 0; // Unless nearest neighbor case is selected (ignored)

        switch (resampleType) {
            case Interpolation.INTERP_NEAREST: // Return a copy of the qsFilter
                return (float[]) qsFilter.clone();
            case Interpolation.INTERP_BILINEAR: // 2 by 2 resample filter
                float[] bilinearKernel = {1.0F / 2.0F};
                return convolveSymmetricKernels(qsParity, resampParity, qsFilter, bilinearKernel);
            case Interpolation.INTERP_BICUBIC: // 4 by 4 resample filter
                float[] bicubicKernel = {9.0F / 16.0F, -1.0F / 16.0F};
                return convolveSymmetricKernels(qsParity, resampParity, qsFilter, bicubicKernel);
            case Interpolation.INTERP_BICUBIC_2: // alternate 4 by 4 resample filter
                float[] bicubic2Kernel = {5.0F / 8.0F, -1.0F / 8.0F};
                return convolveSymmetricKernels(qsParity, resampParity, qsFilter, bicubic2Kernel);
            default:
                throw new IllegalArgumentException(JaiI18N.getString("FilteredSubsample0"));
        }
    } // combineFilters

    /**
     * <code>filterParity</code> -- Returns combined filter/resample parity. This is odd only when we have an odd scale
     * factor or we have nearest neighbor filtering (no resample kernel needed). <code>scaleFactor</code> was validated
     * by the constructor. Possible return values are 0 or 1.
     *
     * @param scaleFactor positive int representing the downsample factor.
     * @param resampleType int representing the interpolation type.
     * @return int representing combined filter parity (0 or 1).
     */
    private static int filterParity(int scaleFactor, int resampleType) {

        // Test scale factor for oddness or nearest neighbor resampling
        if ((scaleFactor % 2 == 1) || (resampleType == Interpolation.INTERP_NEAREST)) return 1;

        // for all other cases we will be convolving an odd filter with an even
        // filter, thus producing an even filter
        return 0;
    } // filterParity

    /**
     * <code>layoutHelper</code> validates input and returns an <code>ImageLayout</code> object.
     *
     * @param source a RenderedImage object.
     * @param interp an Interpolation object.
     * @param scaleX an int downsample factor.
     * @param scaleY an int downsample factor.
     * @param filterSize an int representing the size of the combined filter and resample kernel.
     * @param il an ImageLayout object.
     * @return validated ImageLayout object.
     */
    private static final ImageLayout layoutHelper(
            RenderedImage source, Interpolation interp, int scaleX, int scaleY, int filterSize, ImageLayout il) {

        if (scaleX < 1 || scaleY < 1) {
            throw new IllegalArgumentException(JaiI18N.getString("FilteredSubsample1"));
        }
        if (filterSize < 1) {
            throw new IllegalArgumentException(JaiI18N.getString("FilteredSubsample2"));
        }

        // Set the bounds to the scaled source bounds.
        Rectangle bounds = forwardMapRect(
                source.getMinX(), source.getMinY(), source.getWidth(), source.getHeight(), scaleX, scaleY);

        // If the user has supplied a layout, use it
        ImageLayout layout = (il == null)
                ? new ImageLayout(bounds.x, bounds.y, bounds.width, bounds.height)
                : (ImageLayout) il.clone();

        // Override dimensions if user passed a hint
        if (il != null) {
            layout.setWidth(bounds.width);
            layout.setHeight(bounds.height);
            layout.setMinX(bounds.x);
            layout.setMinY(bounds.y);
        }

        return layout;
    } // setLayout

    /**
     * <code>FilteredSubsampleOpImage</code> constructs an OpImage representing filtered integral subsampling. The scale
     * factors represent the ratio of source to destination dimensions.
     *
     * @param source a RenderedImage.
     * @param extender a BorderExtender, or null.
     * @param config a Map object possibly holding tile cache information
     * @param layout an ImageLayout optionally containing the tile grid layout, SampleModel, and ColorModel, or null.
     * @param interp a Interpolation object to use for resampling.
     * @param scaleX downsample factor along x axis.
     * @param scaleY downsample factor along y axis.
     * @param qsFilter symmetric filter coefficients (partial kernel).
     * @throws IllegalArgumentException if the interp type is not one of: INTERP_NEAREST, INTERP_BILINEAR,
     *     INTERP_BICUBIC, or INTERP_BICUBIC_2
     */
    public FilteredSubsampleOpImage(
            RenderedImage source,
            BorderExtender extender,
            Map config,
            ImageLayout layout,
            int scaleX,
            int scaleY,
            float[] qsFilter,
            Interpolation interp) {

        // Propagate to GeometricOpImage constructor
        super(
                vectorize(source),
                layoutHelper(source, interp, scaleX, scaleY, qsFilter.length, layout),
                config, // Map object
                true, // cobbleSources,
                extender, // extender
                interp, // Interpolation object
                null);

        int resampleType;

        // Determine the interpolation type, if not supported throw exception
        if (interp instanceof InterpolationNearest) {
            resampleType = Interpolation.INTERP_NEAREST;
        } else if (interp instanceof InterpolationBilinear) {
            resampleType = Interpolation.INTERP_BILINEAR;
        } else if (interp instanceof InterpolationBicubic) {
            resampleType = Interpolation.INTERP_BICUBIC;
        } else if (interp instanceof InterpolationBicubic2) {
            resampleType = Interpolation.INTERP_BICUBIC_2;
        } else {
            throw new IllegalArgumentException(JaiI18N.getString("FilteredSubsample3"));
        }

        // Construct combined anti-alias and resample kernels.
        this.hParity = filterParity(scaleX, resampleType);
        this.vParity = filterParity(scaleY, resampleType);
        this.hKernel = combineFilters(scaleX, resampleType, qsFilter);
        this.vKernel = combineFilters(scaleY, resampleType, qsFilter);

        this.scaleX = scaleX;
        this.scaleY = scaleY;
    } // FilteredSubsampleOpImage

    /**
     * Computes the source point corresponding to the supplied point.
     *
     * @param destPt the position in destination image coordinates to map to source image coordinates.
     * @return a <code>Point2D</code> of the same class as <code>destPt</code>.
     * @throws IllegalArgumentException if <code>destPt</code> is <code>null</code>.
     * @since JAI 1.1.2
     */
    public Point2D mapDestPoint(Point2D destPt) {
        if (destPt == null) {
            throw new IllegalArgumentException(JaiI18N.getString("Generic0"));
        }

        Point2D pt = (Point2D) destPt.clone();

        pt.setLocation(destPt.getX() * scaleX, destPt.getY() * scaleY);

        return pt;
    }

    /**
     * Computes the destination point corresponding to the supplied point.
     *
     * @param sourcePt the position in source image coordinates to map to destination image coordinates.
     * @return a <code>Point2D</code> of the same class as <code>sourcePt</code>.
     * @throws IllegalArgumentException if <code>sourcePt</code> is <code>null</code>.
     * @since JAI 1.1.2
     */
    public Point2D mapSourcePoint(Point2D sourcePt) {
        if (sourcePt == null) {
            throw new IllegalArgumentException(JaiI18N.getString("Generic0"));
        }

        Point2D pt = (Point2D) sourcePt.clone();

        pt.setLocation(sourcePt.getX() / scaleX, sourcePt.getY() / scaleY);

        return pt;
    }

    /**
     * Returns a conservative estimate of the destination region that can potentially be affected by the pixels of a
     * rectangle of a given source.
     *
     * @param sourceRect The <code>Rectangle</code> in source coordinates.
     * @param sourceIndex The index of the source image.
     * @return a <code>Rectangle</code> indicating the potentially affected destination region, or <code>null</code> if
     *     the region is unknown.
     * @throws IllegalArgumentException if <code>sourceIndex</code> is negative or greater than the index of the last
     *     source.
     * @throws NullPointerException if <code>sourceRect</code> is <code>null</code>.
     */
    public Rectangle mapSourceRect(Rectangle sourceRect, int sourceIndex) {
        if (sourceIndex != 0) { // this image only has one source
            throw new IllegalArgumentException(JaiI18N.getString("FilteredSubsample4"));
        }

        int xOffset = sourceRect.x + hKernel.length - hParity - scaleX / 2;
        int yOffset = sourceRect.y + vKernel.length - vParity - scaleY / 2;
        int rectWidth = sourceRect.width - 2 * hKernel.length + hParity + 1;
        int rectHeight = sourceRect.height - 2 * vKernel.length + vParity + 1;
        return forwardMapRect(xOffset, yOffset, rectWidth, rectHeight, scaleX, scaleY);
    } // mapSourceRect

    /**
     * Forward map a source Rectangle into destination space.
     *
     * @param x source frame coordinate.
     * @param y source frame coordinate.
     * @param w source frame width.
     * @param h source frame height.
     * @param scaleX downsample factor.
     * @param scaleY downsample factor.
     * @return a <code>Rectangle</code> indicating the destination region.
     */
    private static final Rectangle forwardMapRect(int x, int y, int w, int h, int scaleX, int scaleY) {
        float sx = 1.0F / scaleX;
        float sy = 1.0F / scaleY;

        x = Math.round(x * sx);
        y = Math.round(y * sy);

        return new Rectangle(x, y, Math.round((x + w) * sx) - x, Math.round((y + h) * sy) - y);
    } // forwardMapRect

    /**
     * Forward map a source Rectangle into destination space. Required by abstract GeometricOpImage
     *
     * @param srcRect a source Rectangle.
     * @param srcIndex int source index (0 for this operator)
     * @return a <code>Rectangle</code> indicating the destination region.
     */
    protected final Rectangle forwardMapRect(Rectangle srcRect, int srcIndex) {
        int x = srcRect.x;
        int y = srcRect.y;
        int w = srcRect.width;
        int h = srcRect.height;
        float sx = 1.0F / scaleX;
        float sy = 1.0F / scaleY;

        x = Math.round(x * sx);
        y = Math.round(y * sy);

        return new Rectangle(x, y, Math.round((x + w) * sx) - x, Math.round((y + h) * sy) - y);
    } // forwardMapRect

    /**
     * Backward map a destination Rectangle into source space.
     *
     * @param destRect a destination Rectangle.
     * @param srcIndex int source index (0 for this operator)
     * @return a <code>Rectangle</code> indicating the source region.
     */
    protected final Rectangle backwardMapRect(Rectangle destRect, int srcIncex) {
        int x = destRect.x;
        int y = destRect.y;
        int w = destRect.width;
        int h = destRect.height;

        return new Rectangle(x * scaleX, y * scaleY, (x + w) * scaleX - x, (y + h) * scaleY - y);
    } // backwardMapRect

    /**
     * Returns a conservative estimate of the region of a specified source that is required in order to compute the
     * pixels of a given destination rectangle.
     *
     * @param destRect The <code>Rectangle</code> in destination coordinates.
     * @param sourceIndex The index of the source image.
     * @return a <code>Rectangle</code> indicating the required source region.
     * @throws IllegalArgumentException if <code>sourceIndex</code> is negative or greater than the index of the last
     *     source.
     * @throws NullPointerException if <code>destRect</code> is <code>null</code>.
     */
    public Rectangle mapDestRect(Rectangle destRect, int sourceIndex) {
        if (sourceIndex != 0) { // this image only has one source
            throw new IllegalArgumentException(JaiI18N.getString("FilteredSubsample4"));
        }
        int xOffset = destRect.x * scaleX - hKernel.length + hParity + scaleX / 2;
        int yOffset = destRect.y * scaleY - vKernel.length + vParity + scaleY / 2;
        int rectWidth = destRect.width * scaleX + 2 * hKernel.length - hParity - 1;
        int rectHeight = destRect.height * scaleY + 2 * vKernel.length - vParity - 1;
        return new Rectangle(xOffset, yOffset, rectWidth, rectHeight);
    } // mapDestRect

    /**
     * Performs a combined subsample/filter operation on a specified rectangle. The sources are cobbled.
     *
     * @param sources an array of source Rasters, guaranteed to provide all necessary source data for computing the
     *     output.
     * @param dest a WritableRaster containing the area to be computed.
     * @param destRect the rectangle within dest to be processed.
     */
    public void computeRect(Raster[] sources, WritableRaster dest, Rectangle destRect) {

        // Get RasterAccessor tags (initialized in OpImage superclass).
        RasterFormatTag[] formatTags = getFormatTags();

        // Get destination accessor.
        RasterAccessor dst = new RasterAccessor(dest, destRect, formatTags[1], getColorModel());

        // Get source accessor.
        RasterAccessor src = new RasterAccessor(
                sources[0],
                mapDestRect(destRect, 0),
                formatTags[0],
                getSourceImage(0).getColorModel());

        switch (dst.getDataType()) {
            case DataBuffer.TYPE_BYTE:
                computeRectByte(src, dst);
                break;
            case DataBuffer.TYPE_USHORT:
                computeRectUShort(src, dst);
                break;
            case DataBuffer.TYPE_SHORT:
                computeRectShort(src, dst);
                break;
            case DataBuffer.TYPE_INT:
                computeRectInt(src, dst);
                break;
            case DataBuffer.TYPE_FLOAT:
                computeRectFloat(src, dst);
                break;
            case DataBuffer.TYPE_DOUBLE:
                computeRectDouble(src, dst);
                break;
            default:
                throw new IllegalArgumentException(JaiI18N.getString("FilteredSubsample5"));
        }

        // If the RasterAccessor set up a temporary write buffer for the
        // operator, tell it to copy that data to the destination Raster.
        if (dst.isDataCopy()) {
            dst.clampDataArrays();
            dst.copyDataToRaster();
        }
    } // computeRect

    /**
     * <code>computeRectByte</code> filter subsamples byte pixel data.
     *
     * @param src RasterAccessor for source image.
     * @param dst RasterAccessor for output image.
     */
    protected void computeRectByte(RasterAccessor src, RasterAccessor dst) {

        // Get dimensions.
        int dwidth = dst.getWidth();
        int dheight = dst.getHeight();
        int dnumBands = dst.getNumBands();

        // Get destination data array references and strides.
        byte dstDataArrays[][] = dst.getByteDataArrays();
        int dstBandOffsets[] = dst.getBandOffsets();
        int dstPixelStride = dst.getPixelStride();
        int dstScanlineStride = dst.getScanlineStride();

        // Get source data array references and strides.
        byte srcDataArrays[][] = src.getByteDataArrays();
        int srcBandOffsets[] = src.getBandOffsets();
        int srcPixelStride = src.getPixelStride();
        int srcScanlineStride = src.getScanlineStride();

        // Compute reused numbers
        int kernelNx = 2 * hKernel.length - hParity;
        int kernelNy = 2 * vKernel.length - vParity;
        int stepDown = (kernelNy - 1) * srcScanlineStride;
        int stepRight = (kernelNx - 1) * srcPixelStride;

        int upLeft, upRight, dnLeft, dnRight;
        float kk;
        float vCtr = vKernel[0];
        float hCtr = hKernel[0];
        int kInd, xLeft, xRight, xUp, xDown;

        for (int band = 0; band < dnumBands; band++) {
            byte dstData[] = dstDataArrays[band];
            byte srcData[] = srcDataArrays[band];
            int srcScanlineOffset = srcBandOffsets[band];
            int dstScanlineOffset = dstBandOffsets[band];

            // Step over source raster coordinates
            for (int ySrc = 0; ySrc < scaleY * dheight; ySrc += scaleY) {
                int dInd = dstScanlineOffset;
                for (int xSrc = 0; xSrc < scaleX * dwidth; xSrc += scaleX) {

                    int upLeft0 = xSrc * srcPixelStride + ySrc * srcScanlineStride + srcScanlineOffset;
                    int upRight0 = upLeft0 + stepRight;
                    int dnLeft0 = upLeft0 + stepDown;
                    int dnRight0 = upRight0 + stepDown;

                    // Exploit 4-fold symmetry
                    float sum = 0;

                    // Make the rectangle squeeze in the vertical direction
                    for (int iy = vKernel.length - 1; iy > vParity - 1; iy--) {
                        upLeft = upLeft0;
                        upRight = upRight0;
                        dnLeft = dnLeft0;
                        dnRight = dnRight0;

                        // Make the rectangle squeeze in the horizontal direction
                        for (int ix = hKernel.length - 1; ix > hParity - 1; ix--) {
                            kk = hKernel[ix] * vKernel[iy];
                            sum += kk
                                    * ((int) (srcData[upLeft] & 0xff)
                                            + (int) (srcData[upRight] & 0xff)
                                            + (int) (srcData[dnLeft] & 0xff)
                                            + (int) (srcData[dnRight] & 0xff));
                            upLeft += srcPixelStride;
                            upRight -= srcPixelStride;
                            dnLeft += srcPixelStride;
                            dnRight -= srcPixelStride;
                        } // ix
                        upLeft0 += srcScanlineStride; // down a row
                        upRight0 += srcScanlineStride;
                        dnLeft0 -= srcScanlineStride; // up a row
                        dnRight0 -= srcScanlineStride;
                    } // iy

                    // Compute the remaining 2-Fold symmetry portions (across and down as needed)

                    // Loop down the center (hParity is odd)
                    if (hParity == 1) {
                        xUp = (xSrc + hKernel.length - 1) * srcPixelStride
                                + ySrc * srcScanlineStride
                                + srcScanlineOffset;
                        xDown = xUp + stepDown;
                        kInd = vKernel.length - 1;
                        while (xUp < xDown) {
                            kk = hCtr * vKernel[kInd--];
                            sum += kk * ((int) (srcData[xUp] & 0xff) + (int) (srcData[xDown] & 0xff));
                            xUp += srcScanlineStride;
                            xDown -= srcScanlineStride;
                        }
                    } // hParity

                    // Loop across the center (vParity is odd), pick up the center if hParity was odd
                    if (vParity == 1) {
                        xLeft = xSrc * srcPixelStride
                                + (ySrc + vKernel.length - 1) * srcScanlineStride
                                + srcScanlineOffset;
                        xRight = xLeft + stepRight;
                        kInd = hKernel.length - 1;
                        while (xLeft < xRight) {
                            kk = vCtr * hKernel[kInd--];
                            sum += kk * ((int) (srcData[xLeft] & 0xff) + (int) (srcData[xRight] & 0xff));
                            xLeft += srcPixelStride;
                            xRight -= srcPixelStride;
                        } // while xLeft

                        // Grab the center pixel if hParity was odd also
                        if (hParity == 1) sum += vCtr * hCtr * (int) (srcData[xLeft] & 0xff);
                    } // if vParity

                    // Convert the sum to an output pixel
                    if (sum < 0.0) sum = 0;
                    if (sum > 255.0) sum = 255;

                    dstData[dInd] = (byte) (sum + 0.5);

                    dInd += dstPixelStride;
                } // for xSrc

                dstScanlineOffset += dstScanlineStride;
            } // for ySrc
        } // for band

        return;
    } // computeRectByte

    /**
     * <code>computeRectUShort</code> filter subsamples unsigned short pixel data.
     *
     * @param src RasterAccessor for source image.
     * @param dst RasterAccessor for output image.
     */
    protected void computeRectUShort(RasterAccessor src, RasterAccessor dst) {

        // Get dimensions.
        int dwidth = dst.getWidth();
        int dheight = dst.getHeight();
        int dnumBands = dst.getNumBands();

        // Get destination data array references and strides.
        short dstDataArrays[][] = dst.getShortDataArrays();
        int dstBandOffsets[] = dst.getBandOffsets();
        int dstPixelStride = dst.getPixelStride();
        int dstScanlineStride = dst.getScanlineStride();

        // Get source data array references and strides.
        short srcDataArrays[][] = src.getShortDataArrays();
        int srcBandOffsets[] = src.getBandOffsets();
        int srcPixelStride = src.getPixelStride();
        int srcScanlineStride = src.getScanlineStride();

        // Compute reused numbers
        int kernelNx = 2 * hKernel.length - hParity;
        int kernelNy = 2 * vKernel.length - vParity;
        int stepDown = (kernelNy - 1) * srcScanlineStride;
        int stepRight = (kernelNx - 1) * srcPixelStride;

        int upLeft, upRight, dnLeft, dnRight;
        float kk;
        float vCtr = vKernel[0];
        float hCtr = hKernel[0];
        int kInd, xLeft, xRight, xUp, xDown;

        for (int band = 0; band < dnumBands; band++) {
            short dstData[] = dstDataArrays[band];
            short srcData[] = srcDataArrays[band];
            int srcScanlineOffset = srcBandOffsets[band];
            int dstScanlineOffset = dstBandOffsets[band];

            // Step over source raster coordinates
            for (int ySrc = 0; ySrc < scaleY * dheight; ySrc += scaleY) {
                int dInd = dstScanlineOffset;
                for (int xSrc = 0; xSrc < scaleX * dwidth; xSrc += scaleX) {

                    int upLeft0 = xSrc * srcPixelStride + ySrc * srcScanlineStride + srcScanlineOffset;
                    int upRight0 = upLeft0 + stepRight;
                    int dnLeft0 = upLeft0 + stepDown;
                    int dnRight0 = upRight0 + stepDown;

                    // Exploit 4-fold symmetry
                    float sum = 0;

                    // Make the rectangle squeeze in the vertical direction
                    for (int iy = vKernel.length - 1; iy > vParity - 1; iy--) {
                        upLeft = upLeft0;
                        upRight = upRight0;
                        dnLeft = dnLeft0;
                        dnRight = dnRight0;

                        // Make the rectangle squeeze in the horizontal direction
                        for (int ix = hKernel.length - 1; ix > hParity - 1; ix--) {
                            kk = hKernel[ix] * vKernel[iy];
                            sum += kk
                                    * ((int) (srcData[upLeft] & 0xffff)
                                            + (int) (srcData[upRight] & 0xffff)
                                            + (int) (srcData[dnLeft] & 0xffff)
                                            + (int) (srcData[dnRight] & 0xffff));
                            upLeft += srcPixelStride;
                            upRight -= srcPixelStride;
                            dnLeft += srcPixelStride;
                            dnRight -= srcPixelStride;
                        } // ix
                        upLeft0 += srcScanlineStride; // down a row
                        upRight0 += srcScanlineStride;
                        dnLeft0 -= srcScanlineStride; // up a row
                        dnRight0 -= srcScanlineStride;
                    } // iy

                    // Compute the remaining 2-Fold symmetry portions
                    // (across and down as needed)

                    // Loop down the center (hParity is odd)
                    if (hParity == 1) {
                        xUp = (xSrc + hKernel.length - 1) * srcPixelStride
                                + ySrc * srcScanlineStride
                                + srcScanlineOffset;
                        xDown = xUp + stepDown;
                        kInd = vKernel.length - 1;
                        while (xUp < xDown) {
                            kk = hCtr * vKernel[kInd--];
                            sum += kk * ((int) (srcData[xUp] & 0xffff) + (int) (srcData[xDown] & 0xffff));
                            xUp += srcScanlineStride;
                            xDown -= srcScanlineStride;
                        }
                    } // hParity

                    // Loop across the center (vParity is odd), pick up the center if hParity was odd
                    if (vParity == 1) {
                        xLeft = xSrc * srcPixelStride
                                + (ySrc + vKernel.length - 1) * srcScanlineStride
                                + srcScanlineOffset;
                        xRight = xLeft + stepRight;
                        kInd = hKernel.length - 1;
                        while (xLeft < xRight) {
                            kk = vCtr * hKernel[kInd--];
                            sum += kk * ((int) (srcData[xLeft] & 0xffff) + (int) (srcData[xRight] & 0xffff));
                            xLeft += srcPixelStride;
                            xRight -= srcPixelStride;
                        } // while xLeft

                        // Grab the center pixel if hParity was odd also
                        if (hParity == 1) sum += vCtr * hCtr * (int) (srcData[xLeft] & 0xffff);
                    } // if vParity
                    int val = (int) (sum + 0.5);
                    dstData[dInd] = (short) (val > 0xffff ? 0xffff : (val < 0 ? 0 : val));

                    dInd += dstPixelStride;
                } // for xSrc

                dstScanlineOffset += dstScanlineStride;
            } // for ySrc
        } // for band

        return;
    } // computeRectUShort

    /**
     * <code>computeRectShort</code> filter subsamples short pixel data.
     *
     * @param src RasterAccessor for source image.
     * @param dst RasterAccessor for output image.
     */
    protected void computeRectShort(RasterAccessor src, RasterAccessor dst) {

        // Get dimensions.
        int dwidth = dst.getWidth();
        int dheight = dst.getHeight();
        int dnumBands = dst.getNumBands();

        // Get destination data array references and strides.
        short dstDataArrays[][] = dst.getShortDataArrays();
        int dstBandOffsets[] = dst.getBandOffsets();
        int dstPixelStride = dst.getPixelStride();
        int dstScanlineStride = dst.getScanlineStride();

        // Get source data array references and strides.
        short srcDataArrays[][] = src.getShortDataArrays();
        int srcBandOffsets[] = src.getBandOffsets();
        int srcPixelStride = src.getPixelStride();
        int srcScanlineStride = src.getScanlineStride();

        // Compute reused numbers
        int kernelNx = 2 * hKernel.length - hParity;
        int kernelNy = 2 * vKernel.length - vParity;
        int stepDown = (kernelNy - 1) * srcScanlineStride;
        int stepRight = (kernelNx - 1) * srcPixelStride;

        int upLeft, upRight, dnLeft, dnRight;
        float kk;
        float vCtr = vKernel[0];
        float hCtr = hKernel[0];
        int kInd, xLeft, xRight, xUp, xDown;

        for (int band = 0; band < dnumBands; band++) {
            short dstData[] = dstDataArrays[band];
            short srcData[] = srcDataArrays[band];
            int srcScanlineOffset = srcBandOffsets[band];
            int dstScanlineOffset = dstBandOffsets[band];

            // Step over source raster coordinates
            for (int ySrc = 0; ySrc < scaleY * dheight; ySrc += scaleY) {
                int dInd = dstScanlineOffset;
                for (int xSrc = 0; xSrc < scaleX * dwidth; xSrc += scaleX) {

                    int upLeft0 = xSrc * srcPixelStride + ySrc * srcScanlineStride + srcScanlineOffset;
                    int upRight0 = upLeft0 + stepRight;
                    int dnLeft0 = upLeft0 + stepDown;
                    int dnRight0 = upRight0 + stepDown;

                    // Exploit 4-fold symmetry
                    float sum = 0;

                    // Make the rectangle squeeze in the vertical direction
                    for (int iy = vKernel.length - 1; iy > vParity - 1; iy--) {
                        upLeft = upLeft0;
                        upRight = upRight0;
                        dnLeft = dnLeft0;
                        dnRight = dnRight0;

                        // Make the rectangle squeeze in the horizontal direction
                        for (int ix = hKernel.length - 1; ix > hParity - 1; ix--) {
                            kk = hKernel[ix] * vKernel[iy];
                            sum += kk
                                    * ((int) (srcData[upLeft])
                                            + (int) (srcData[upRight])
                                            + (int) (srcData[dnLeft])
                                            + (int) (srcData[dnRight]));
                            upLeft += srcPixelStride;
                            upRight -= srcPixelStride;
                            dnLeft += srcPixelStride;
                            dnRight -= srcPixelStride;
                        } // ix
                        upLeft0 += srcScanlineStride; // down a row
                        upRight0 += srcScanlineStride;
                        dnLeft0 -= srcScanlineStride; // up a row
                        dnRight0 -= srcScanlineStride;
                    } // iy

                    // Compute the remaining 2-Fold symmetry portions
                    // (across and down as needed)

                    // Loop down the center (hParity is odd)
                    if (hParity == 1) {
                        xUp = (xSrc + hKernel.length - 1) * srcPixelStride
                                + ySrc * srcScanlineStride
                                + srcScanlineOffset;
                        xDown = xUp + stepDown;
                        kInd = vKernel.length - 1;
                        while (xUp < xDown) {
                            kk = hCtr * vKernel[kInd--];
                            sum += kk * ((int) (srcData[xUp]) + (int) (srcData[xDown]));
                            xUp += srcScanlineStride;
                            xDown -= srcScanlineStride;
                        }
                    } // hParity

                    // Loop across the center (vParity is odd), pick up the center if hParity was odd
                    if (vParity == 1) {
                        xLeft = xSrc * srcPixelStride
                                + (ySrc + vKernel.length - 1) * srcScanlineStride
                                + srcScanlineOffset;
                        xRight = xLeft + stepRight;
                        kInd = hKernel.length - 1;
                        while (xLeft < xRight) {
                            kk = vCtr * hKernel[kInd--];
                            sum += kk * ((int) (srcData[xLeft]) + (int) (srcData[xRight]));
                            xLeft += srcPixelStride;
                            xRight -= srcPixelStride;
                        } // while xLeft

                        // Grab the center pixel if hParity was odd also
                        if (hParity == 1) sum += vCtr * hCtr * (int) (srcData[xLeft]);
                    } // if vParity

                    dstData[dInd] = ImageUtil.clampShort((int) (sum + 0.5));
                    dInd += dstPixelStride;
                } // for xSrc

                dstScanlineOffset += dstScanlineStride;
            } // for ySrc
        } // for band

        return;
    } // computeRectShort

    /**
     * <code>computeRectInt</code> filter subsamples int pixel data.
     *
     * @param src RasterAccessor for source image.
     * @param dst RasterAccessor for output image.
     */
    protected void computeRectInt(RasterAccessor src, RasterAccessor dst) {

        // Get dimensions.
        int dwidth = dst.getWidth();
        int dheight = dst.getHeight();
        int dnumBands = dst.getNumBands();

        // Get destination data array references and strides.
        int dstDataArrays[][] = dst.getIntDataArrays();
        int dstBandOffsets[] = dst.getBandOffsets();
        int dstPixelStride = dst.getPixelStride();
        int dstScanlineStride = dst.getScanlineStride();

        // Get source data array references and strides.
        int srcDataArrays[][] = src.getIntDataArrays();
        int srcBandOffsets[] = src.getBandOffsets();
        int srcPixelStride = src.getPixelStride();
        int srcScanlineStride = src.getScanlineStride();

        // Compute reused numbers
        int kernelNx = 2 * hKernel.length - hParity;
        int kernelNy = 2 * vKernel.length - vParity;
        int stepDown = (kernelNy - 1) * srcScanlineStride;
        int stepRight = (kernelNx - 1) * srcPixelStride;

        int upLeft, upRight, dnLeft, dnRight;
        double kk;
        double vCtr = (double) vKernel[0];
        double hCtr = (double) hKernel[0];
        int kInd, xLeft, xRight, xUp, xDown;

        for (int band = 0; band < dnumBands; band++) {
            int dstData[] = dstDataArrays[band];
            int srcData[] = srcDataArrays[band];
            int srcScanlineOffset = srcBandOffsets[band];
            int dstScanlineOffset = dstBandOffsets[band];

            // Step over source raster coordinates
            for (int ySrc = 0; ySrc < scaleY * dheight; ySrc += scaleY) {
                int dInd = dstScanlineOffset;
                for (int xSrc = 0; xSrc < scaleX * dwidth; xSrc += scaleX) {

                    int upLeft0 = xSrc * srcPixelStride + ySrc * srcScanlineStride + srcScanlineOffset;
                    int upRight0 = upLeft0 + stepRight;
                    int dnLeft0 = upLeft0 + stepDown;
                    int dnRight0 = upRight0 + stepDown;

                    // Exploit 4-fold symmetry
                    double sum = 0;

                    // Make the rectangle squeeze in the vertical direction
                    for (int iy = vKernel.length - 1; iy > vParity - 1; iy--) {
                        upLeft = upLeft0;
                        upRight = upRight0;
                        dnLeft = dnLeft0;
                        dnRight = dnRight0;

                        // Make the rectangle squeeze in the horizontal direction
                        for (int ix = hKernel.length - 1; ix > hParity - 1; ix--) {
                            kk = hKernel[ix] * vKernel[iy];
                            sum += kk
                                    * ((long) srcData[upLeft]
                                            + (long) srcData[upRight]
                                            + (long) srcData[dnLeft]
                                            + (long) srcData[dnRight]);
                            upLeft += srcPixelStride;
                            upRight -= srcPixelStride;
                            dnLeft += srcPixelStride;
                            dnRight -= srcPixelStride;
                        } // ix
                        upLeft0 += srcScanlineStride; // down a row
                        upRight0 += srcScanlineStride;
                        dnLeft0 -= srcScanlineStride; // up a row
                        dnRight0 -= srcScanlineStride;
                    } // iy

                    // Compute the remaining 2-Fold symmetry portions
                    // (across and down as needed)

                    // Loop down the center (hParity is odd)
                    if (hParity == 1) {
                        xUp = (xSrc + hKernel.length - 1) * srcPixelStride
                                + ySrc * srcScanlineStride
                                + srcScanlineOffset;
                        xDown = xUp + stepDown;
                        kInd = vKernel.length - 1;
                        while (xUp < xDown) {
                            kk = hCtr * vKernel[kInd--];
                            sum += kk * ((long) srcData[xUp] + (long) srcData[xDown]);
                            xUp += srcScanlineStride;
                            xDown -= srcScanlineStride;
                        }
                    } // hParity

                    // Loop across the center (vParity is odd), pick up the center if hParity was odd
                    if (vParity == 1) {
                        xLeft = xSrc * srcPixelStride
                                + (ySrc + vKernel.length - 1) * srcScanlineStride
                                + srcScanlineOffset;
                        xRight = xLeft + stepRight;
                        kInd = hKernel.length - 1;
                        while (xLeft < xRight) {
                            kk = vCtr * hKernel[kInd--];
                            sum += kk * ((long) (srcData[xLeft]) + (long) (srcData[xRight]));
                            xLeft += srcPixelStride;
                            xRight -= srcPixelStride;
                        } // while xLeft

                        // Grab the center pixel if hParity was odd also
                        if (hParity == 1) sum += vCtr * hCtr * ((long) srcData[xLeft]);
                    } // if vParity

                    dstData[dInd] = ImageUtil.clampInt((int) (sum + 0.5));

                    dInd += dstPixelStride;
                } // for xSrc

                dstScanlineOffset += dstScanlineStride;
            } // for ySrc
        } // for band

        return;
    } // computeRectInt

    /**
     * <code>computeRectFloat</code> filter subsamples float pixel data.
     *
     * @param src RasterAccessor for source image.
     * @param dst RasterAccessor for output image.
     */
    protected void computeRectFloat(RasterAccessor src, RasterAccessor dst) {

        // Get dimensions.
        int dwidth = dst.getWidth();
        int dheight = dst.getHeight();
        int dnumBands = dst.getNumBands();

        // Get destination data array references and strides.
        float dstDataArrays[][] = dst.getFloatDataArrays();
        int dstBandOffsets[] = dst.getBandOffsets();
        int dstPixelStride = dst.getPixelStride();
        int dstScanlineStride = dst.getScanlineStride();

        // Get source data array references and strides.
        float srcDataArrays[][] = src.getFloatDataArrays();
        int srcBandOffsets[] = src.getBandOffsets();
        int srcPixelStride = src.getPixelStride();
        int srcScanlineStride = src.getScanlineStride();

        // Compute reused numbers
        int kernelNx = 2 * hKernel.length - hParity;
        int kernelNy = 2 * vKernel.length - vParity;
        int stepDown = (kernelNy - 1) * srcScanlineStride;
        int stepRight = (kernelNx - 1) * srcPixelStride;

        int upLeft, upRight, dnLeft, dnRight;
        double kk;
        double vCtr = (double) vKernel[0];
        double hCtr = (double) hKernel[0];
        int kInd, xLeft, xRight, xUp, xDown;

        for (int band = 0; band < dnumBands; band++) {
            float dstData[] = dstDataArrays[band];
            float srcData[] = srcDataArrays[band];
            int srcScanlineOffset = srcBandOffsets[band];
            int dstScanlineOffset = dstBandOffsets[band];

            // Step over source raster coordinates
            for (int ySrc = 0; ySrc < scaleY * dheight; ySrc += scaleY) {
                int dInd = dstScanlineOffset;
                for (int xSrc = 0; xSrc < scaleX * dwidth; xSrc += scaleX) {

                    int upLeft0 = xSrc * srcPixelStride + ySrc * srcScanlineStride + srcScanlineOffset;
                    int upRight0 = upLeft0 + stepRight;
                    int dnLeft0 = upLeft0 + stepDown;
                    int dnRight0 = upRight0 + stepDown;

                    // Exploit 4-fold symmetry
                    double sum = 0;

                    // Make the rectangle squeeze in the vertical direction
                    for (int iy = vKernel.length - 1; iy > vParity - 1; iy--) {
                        upLeft = upLeft0;
                        upRight = upRight0;
                        dnLeft = dnLeft0;
                        dnRight = dnRight0;

                        // Make the rectangle squeeze in the horizontal direction
                        for (int ix = hKernel.length - 1; ix > hParity - 1; ix--) {
                            kk = hKernel[ix] * vKernel[iy];
                            sum += kk
                                    * ((double) srcData[upLeft]
                                            + (double) srcData[upRight]
                                            + (double) srcData[dnLeft]
                                            + (double) srcData[dnRight]);
                            upLeft += srcPixelStride;
                            upRight -= srcPixelStride;
                            dnLeft += srcPixelStride;
                            dnRight -= srcPixelStride;
                        } // ix
                        upLeft0 += srcScanlineStride; // down a row
                        upRight0 += srcScanlineStride;
                        dnLeft0 -= srcScanlineStride; // up a row
                        dnRight0 -= srcScanlineStride;
                    } // iy

                    // Compute the remaining 2-Fold symmetry portions
                    // (across and down as needed)

                    // Loop down the center (hParity is odd)
                    if (hParity == 1) {
                        xUp = (xSrc + hKernel.length - 1) * srcPixelStride
                                + ySrc * srcScanlineStride
                                + srcScanlineOffset;
                        xDown = xUp + stepDown;
                        kInd = vKernel.length - 1;
                        while (xUp < xDown) {
                            kk = hCtr * vKernel[kInd--];
                            sum += kk * ((double) srcData[xUp] + (double) srcData[xDown]);
                            xUp += srcScanlineStride;
                            xDown -= srcScanlineStride;
                        }
                    } // hParity

                    // Loop across the center (vParity is odd), pick up the center if hParity was odd
                    if (vParity == 1) {
                        xLeft = xSrc * srcPixelStride
                                + (ySrc + vKernel.length - 1) * srcScanlineStride
                                + srcScanlineOffset;
                        xRight = xLeft + stepRight;
                        kInd = hKernel.length - 1;
                        while (xLeft < xRight) {
                            kk = vCtr * hKernel[kInd--];
                            sum += kk * ((double) (srcData[xLeft]) + (double) (srcData[xRight]));
                            xLeft += srcPixelStride;
                            xRight -= srcPixelStride;
                        } // while xLeft

                        // Grab the center pixel if hParity was odd also
                        if (hParity == 1) sum += vCtr * hCtr * ((double) srcData[xLeft]);
                    } // if vParity

                    dstData[dInd] = ImageUtil.clampFloat(sum);

                    dInd += dstPixelStride;
                } // for xSrc

                dstScanlineOffset += dstScanlineStride;
            } // for ySrc
        } // for band

        return;
    } // computeRectFloat

    /**
     * <code>computeRectDouble</code> filter subsamples double pixel data.
     *
     * @param src RasterAccessor for source image.
     * @param dst RasterAccessor for output image.
     */
    protected void computeRectDouble(RasterAccessor src, RasterAccessor dst) {

        // Get dimensions.
        int dwidth = dst.getWidth();
        int dheight = dst.getHeight();
        int dnumBands = dst.getNumBands();

        // Get destination data array references and strides.
        double dstDataArrays[][] = dst.getDoubleDataArrays();
        int dstBandOffsets[] = dst.getBandOffsets();
        int dstPixelStride = dst.getPixelStride();
        int dstScanlineStride = dst.getScanlineStride();

        // Get source data array references and strides.
        double srcDataArrays[][] = src.getDoubleDataArrays();
        int srcBandOffsets[] = src.getBandOffsets();
        int srcPixelStride = src.getPixelStride();
        int srcScanlineStride = src.getScanlineStride();

        // Compute reused numbers
        int kernelNx = 2 * hKernel.length - hParity;
        int kernelNy = 2 * vKernel.length - vParity;
        int stepDown = (kernelNy - 1) * srcScanlineStride;
        int stepRight = (kernelNx - 1) * srcPixelStride;

        int upLeft, upRight, dnLeft, dnRight;
        double kk;
        double vCtr = (double) vKernel[0];
        double hCtr = (double) hKernel[0];
        int kInd, xLeft, xRight, xUp, xDown;

        for (int band = 0; band < dnumBands; band++) {
            double dstData[] = dstDataArrays[band];
            double srcData[] = srcDataArrays[band];
            int srcScanlineOffset = srcBandOffsets[band];
            int dstScanlineOffset = dstBandOffsets[band];

            // Step over source raster coordinates
            for (int ySrc = 0; ySrc < scaleY * dheight; ySrc += scaleY) {
                int dInd = dstScanlineOffset;
                for (int xSrc = 0; xSrc < scaleX * dwidth; xSrc += scaleX) {

                    int upLeft0 = xSrc * srcPixelStride + ySrc * srcScanlineStride + srcScanlineOffset;
                    int upRight0 = upLeft0 + stepRight;
                    int dnLeft0 = upLeft0 + stepDown;
                    int dnRight0 = upRight0 + stepDown;

                    // Exploit 4-fold symmetry
                    double sum = 0;

                    // Make the rectangle squeeze in the vertical direction
                    for (int iy = vKernel.length - 1; iy > vParity - 1; iy--) {
                        upLeft = upLeft0;
                        upRight = upRight0;
                        dnLeft = dnLeft0;
                        dnRight = dnRight0;

                        // Make the rectangle squeeze in the horizontal direction
                        for (int ix = hKernel.length - 1; ix > hParity - 1; ix--) {
                            kk = hKernel[ix] * vKernel[iy];
                            sum += kk * (srcData[upLeft] + srcData[upRight] + srcData[dnLeft] + srcData[dnRight]);
                            upLeft += srcPixelStride;
                            upRight -= srcPixelStride;
                            dnLeft += srcPixelStride;
                            dnRight -= srcPixelStride;
                        } // ix
                        upLeft0 += srcScanlineStride; // down a row
                        upRight0 += srcScanlineStride;
                        dnLeft0 -= srcScanlineStride; // up a row
                        dnRight0 -= srcScanlineStride;
                    } // iy

                    // Compute the remaining 2-Fold symmetry portions
                    // (across and down as needed)

                    // Loop down the center (hParity is odd)
                    if (hParity == 1) {
                        xUp = (xSrc + hKernel.length - 1) * srcPixelStride
                                + ySrc * srcScanlineStride
                                + srcScanlineOffset;
                        xDown = xUp + stepDown;
                        kInd = vKernel.length - 1;
                        while (xUp < xDown) {
                            kk = hCtr * vKernel[kInd--];
                            sum += kk * (srcData[xUp] + srcData[xDown]);
                            xUp += srcScanlineStride;
                            xDown -= srcScanlineStride;
                        }
                    } // hParity

                    // Loop across the center (vParity is odd), pick up the center if hParity was odd
                    if (vParity == 1) {
                        xLeft = xSrc * srcPixelStride
                                + (ySrc + vKernel.length - 1) * srcScanlineStride
                                + srcScanlineOffset;
                        xRight = xLeft + stepRight;
                        kInd = hKernel.length - 1;
                        while (xLeft < xRight) {
                            kk = vCtr * hKernel[kInd--];
                            sum += kk * (srcData[xLeft] + srcData[xRight]);
                            xLeft += srcPixelStride;
                            xRight -= srcPixelStride;
                        } // while xLeft

                        // Grab the center pixel if hParity was odd also
                        if (hParity == 1) sum += vCtr * hCtr * srcData[xLeft];
                    } // if vParity

                    dstData[dInd] = sum;

                    dInd += dstPixelStride;
                } // for xSrc

                dstScanlineOffset += dstScanlineStride;
            } // for ySrc
        } // for band

        return;
    } // computeRectDouble
} // class FilteredSubsampleOpImage
